function model_settings = model_ini(dx,dt,mode,FitTo,FittingDataLength,vargin)
% MODEL_INI initialize model. Run once before calling MPC1xxb.m
% 
% model_settings = model_ini(dx,dt,mode,FitTo,FittingDataLength) creates a structure with
% information describing the simulation model. The identification table is
% kept empty.
%
% model_settings = model_ini(dx,dt,mode,FileName) 
% 
% INPUT:
%   FileName: The name of the file that save the identified coefficient
%   dx: spacial discretization interval.
%   dt: sampling time.
%   mode: number of mode used.
%
% OUTPUT:
%   model_settings: a structure with the following fields
%     FitTo  = a string describing the fitting variable
%     FittingTimeRange = time length of fitting data. Unit 'second'.
%     W      = several deposition rates that we run open loop simulations for identification.
%     nu     =
%     sigma2 =
%     K      =
%     Tau    =
%     M2ModeWeighting = Weighting of different modes when calculate root mean
%                       square slope. This should be precalculated using
%                       MODEL_ID.
%     mode   = Number of mode used in decomposed EW equation.
%     dt     = Sampling time.
%     dx     = Spatial discretization interval.
%     L0     = spacial range of EW equation is from [-L0,L0]
%     L      = Number of lattice. L = 2*L0/dx;
%
%% 
if nargin == 3
% Default parameters identified by Gangshi
%     model_settings.W       = [0.10,0.15,0.20,0.50,1.0];
%     model_settings.nu      = [0.0451e-4,1.3915e-05,0.3583e-4,0.1232e-3,0.3939e-3];
%     model_settings.sigma2  = [1.0810e-3,6.2341e-3,2.1983e-2,0.0990,0.3007];
%     model_settings.Rh      = [0.1002,0.1509,0.2021,0.5129,1.0432];
%     model_settings.K	     = [0.99962,0.9955,0.9914,0.9730,0.9532];
%     model_settings.Tau     = [10.0740,7.6978,5.3215,2.1915,1.1525];

% Fitting mean(R2) till t =
%     model_settings.W       = [0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50];
%     model_settings.nu      = [4.8017e-3,3.0642e-3,1.2078e-3,3.3919e-4,2.5736e-4,1.6466e-4,1.2699e-4,7.1933e-5,6.2490e-5];
%     model_settings.sigma2  = [1.7110e-2,1.6759e-2,1.2597e-2,8.3201e-3,8.8603e-3,8.7950e-3,8.9019e-3,8.5683e-3,9.2994e-3];
%     model_settings.Rh      = [1.0020e-1,1.5090e-1,2.0210e-1,2.5390e-1,3.0570e-1,3.5750e-1,4.0930e-1,4.6110e-1,5.1290e-1];
%     model_settings.K       = [0.9996,0.9955,0.9914,0.9883,0.9853,0.9822,0.9791,0.9761,0.9730];
%     model_settings.Tau     = [10.0740,7.6978,5.3215,4.7998,4.2782,3.7565,3.2348,2.7132,2.1915];

% Fitting var(R2) till t = 100s
%     model_settings.FitTo   = 'var(R2)';
%     model_settings.FittingTimeRange = 100;
%     model_settings.W       = [0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50];
%     model_settings.nu      = [7.9731e-4,3.5848e-3,4.6012e-6,2.7075e-5,1.6843e-5,5.7364e-6,-6.9883e-6,-1.8224e-5,1.4013e-5];
%     model_settings.sigma2  = [7.3728e-3,1.5817e-2,4.4131e-3,7.0836e-3,8.564e-3,8.7863e-3,9.1571e-3,9.0314e-3,1.4381e-2];
%     model_settings.Rh      = [0.1002 0.1509 0.2021 0.2539 0.3057 0.3575 0.4093 0.4611 0.5129];
%     model_settings.K       = [0.99962,0.9955 0.9914 0.98833 0.98527 0.9822 0.97913 0.97607 0.973];
%     model_settings.Tau     = [10.074 7.6978 5.3215 4.7998 4.2782 3.7565 3.2348 2.7132 2.1915];

% Fitting var(R2) till t=1000s
model_settings.FitTo   = 'var(R2)';
model_settings.FittingTimeRange = 1000;
model_settings.W       = [0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50];
model_settings.nu      = [3.4941e-4 -9.6277e-7 1.1522e-7 -3.1824e-6 -7.0819e-7 -2.2166e-6 2.6946e-6 -3.374e-7 -5.1585e-7];
model_settings.sigma2  = [6.1298e-3 1.9592e-3 3.4054e-3 3.1876e-3 7.4906e-3 9.6178e-3 1.7936e-2 1.7489e-2 1.9078e-2];
model_settings.Rh      = [0.1002 0.1509 0.2021 0.2539 0.3057 0.3575 0.4093 0.4611 0.5129];
model_settings.K       = [0.99962 0.9955 0.9914 0.98833 0.98527 0.9822 0.97913 0.97607 0.973];
model_settings.Tau     = [10.074 7.6978 5.3215 4.7998 4.2782 3.7565 3.2348 2.7132 2.1915];

% Fitting var(R2) till t=3000s
%     model_settings.W       = [0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50];
%     model_settings.nu      = [2.5874e-7 5.3862e-7 5.8229e-7 4.2839e-7 -5.0617e-7 -4.5717e-7 1.6303e-6 -6.0172e-7 1.6779e-6];
%     model_settings.sigma2  = [0.0013728 0.003022 0.004927 0.0081703 0.0084516 0.012477 0.022967 0.016451 0.031638];
%     model_settings.Rh      = [0.1002 0.1509 0.2021 0.2539 0.3057 0.3575 0.4093 0.4611 0.5129];
%     model_settings.K       = [0.99962 0.9955 0.9914 0.98833 0.98527 0.9822 0.97913 0.97607 0.973];
%     model_settings.Tau     = [10.074 7.6978 5.3215 4.7998 4.2782 3.7565 3.2348
%     2.7132 2.1915];
end

if nargin == 4
    load FileName
    model_settings.FitTo = FitTo;
    model_settings.FittingTimeRange = FittingTimeRange;
    model_settings.W = W;
    model_settings.nu = nu;
    model_settings.sigma2 = sigma2;
    model_settings.Rh = Rh;
    model_settings.K = K;
    model_settings.Tau = Tau;
end

if nargin == 5
    model_settings.FitTo = FitTo;
    model_settings.FittingTimeRange = FittingDataLength;
    model_settings.W = [];
    model_settings.nu = [];
    model_settings.sigma2 = [];
    model_settings.Rh = [];
    model_settings.K = [];
    model_settings.Tau = [];
end
%% 
dx = 2*pi/100;
x = (1:100)*dx-pi;

K_alpha = zeros(mode,1);
K_beta = zeros(mode,1);
for n=1:mode
    h_sin = sin(n*x);
    dh_sin = h_sin;
    dh_sin(1) = h_sin(1)-h_sin(end);
    dh_sin(2:end) = diff(h_sin);
    K_alpha(n) = sum(dh_sin.^2)/(pi*dx^2);

    h_cos = cos(n*x);
    dh_cos = h_cos;
    dh_cos(1) = h_cos(1)-h_cos(end);
    dh_cos(2:end) = diff(h_cos);
    K_beta(n)  = sum(dh_cos.^2)/(pi*dx^2);    
end
model_settings.M2ModeWeighting=[K_alpha,K_beta];

model_settings.L0 = pi;
model_settings.L  = 100;
model_settings.dx = dx;
model_settings.dt = dt;
model_settings.mode = mode;
